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ABSTRACT 

Chaos  describes  a  class  of  motions  of  a  deterministic  system  whose  time 
history  is  sensitive  to  initial  conditions.  Because  of  the  sensitivity  of  initial 
conditions,  the  response  of  a  dynamical  system  may  result  in  instabilities. 
Hence,  a  study  of  nonlinear  response  of  structures  under  the  expected 
frequencies  of  excitation  becomes  important.  Chaotic  behavior,  for  example, 
may  be  found  in  the  vibration  response  of  large  flexible  space  structures 
including  trusses,  booms,  and  radio  antennas.  Methods  of  quantifying  chaos 
have  been  applied  to  flexible  beams  both  analytically  and  experimentally.  This 
research  effort  investigates  the  effects  of  sensors,  strain  gages  and 
accelerometers,  in  studying  chaotic  motions.  A  long  flexible  beam  is  used  to 
model  the  chaotic  behavior,  which  is  also  mathematically  modeled  as  Duffing's 
Equation.  Time  histories  are  recorded  and  analyzed  using  pseudo-phase  space, 
Fourier  spectrums,  Poincare  sections,  Lyapunov  exponents  and  fractal 
correlation  dimensions.   Comparison  of  the  two  sensors  is  also  performed. 
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I.        INTRODUCTION 

A.   REASON   FOR   ANALYSIS 

Many  of  the  satellites  currently  being  launched  have  highly  flexible 
components  such  as  trusses  and  booms.  Long,  highly  flexible  booms  and 
trusses  are  used  for  gravity  gradient  stabilization,  supporting  large  solar  arrays 
and  radio  antennas,  and  ensuring  that  humans  and  equipment  are  adequately 
shielded  from  radioactive  power  generators  in  some  proposed  interplanetary 
spacecraft.  Deploying  these  trusses  and  booms  and  maneuvering  the  satellite 
for  attitude  control  and  station  keeping  result  in  large  displacements  and 
rotations  of  the  tips.  This  leads  to  the  presence  of  geometric  nonlinearities  in 
the  structure.  Material  nonlinearities  are  not  usually  present  since  the  booms 
and  trusses  are  not  deflected  past  their  yield  point  and  thus  the  linear  stress- 
strain  law  applies.  Because  of  the  geometric  nonlinearities,  it  is  appropriate  to 
model  flexible  space  structures  by  nonlinear  equations  of  motions  instead  of 
linear  ones.  Duffing's  Equation  is  one  equation  used  to  adequately  describe  the 
nonlinear  behavior  due  to  the  elastic  effects  of  these  large  structures. 

The  presence  of  nonlinearities  may  result  in  motions  which  are  impossible 
to  predict  after  long  periods  of  time,  yet  completely  deterministic  since  the 
equations  of  motions  are  known.  Slight  changes  in  initial  conditions  produce 
radically  different  outputs.  These  types  of  responses  have  been  characterized 
as  "chaotic."  Chaotic  vibrations  of  spacecraft  structures  could  result  in  control 
problems. 


B.    SCOPE  OF  THESIS 

The  purpose  of  this  theses  is  to  demonstrate  that  chaotic  vibrations  can 
occur  in  highly  flexible  space  structures  using  theoretical  and  experimental 
analyses.  For  the  theoretical  analysis,  a  flexible  sheet  of  aluminum  7  ft.  long,  3 
in.  wide,  and  0.080  in.  thick,  was  modelled  by  Duffing's  Equations.  This  process 
is  described  in  Chapter  III.  For  the  experimental  analysis,  the  same  beam  was 
instrumented  with  strain  gauges  and  accelerometers.  The  beam  was  then 
shaken  by  an  MB  Dynamics  Model  PM500A  Vibration  Exciter  and  the  time 
histories  of  the  strains  and  accelerations  were  recorded.  This  process  is 
discussed  in  Chapter  IV.  Data  from  both  analyses  were  then  analyzed  using 
the  computer  program,  "Chaos,"  developed  by  CDR  Martinus  Sarigul-Klijn, 
USN,  [Ref.  1]. 

In  addition,  the  purpose  of  this  thesis  is  to  investigate  the  sensitivity  of 
sensors  in  detecting  chaotic  vibrations.  The  two  most  commonly  used  types  of 
sensors  are  accelerometers  and  strain  gauges.  Measuring  accelerations  may 
be  a  more  accurate  method  of  recording  the  time  histories  of  chaotic  systems 
since  this  includes  inertia,  damping,  and  stiffness  terms  of  the  mathematical 
equations  used  to  model  the  system.  Strain  gauges  permit  measurement  of  only 
the  stiffness  terms  of  the  mathematical  equations  used  to  model  the  system. 
After  a  discussion  of  the  Science  of  Chaos  (Chapter  II)  the  results  of  the 
theoretical  and  experimental  analyses  will  be  discussed. 


II.      THE  SCIENCE  OF  CHAOS 

A.   BACKGROUND   AND   PREVIOUS   STUDIES 

A  deterministic  system  whose  time  history  has  a  sensitive  dependence  on 
initial  conditions  is  one  definition  of  a  chaotic  system  [Ref.  2:  p.  4].  Even  though 
the  equations  of  motions  which  adequately  describe  a  system  are  well  known, 
small  changes  in  the  initial  conditions  produce  markedly  different  outputs.  For 
example,  although  there  are  many  computer  programs  which  can  accurately 
predict  the  weather  tomorrow,  there  are  no  programs  which  can  accurately 
predict  the  weather  for  a  day  next  year.  Recognizing  that  chaos  occurs  in  the 
nonlinear  deterministic  models  raises  the  hope  of  understanding  the  source  of 
this  random-like  behavior  and  possibly  doing  something  to  change  it. 

Although  Henri  Poincare  observed  chaotic  behavior  in  the  trajectories  of 
some  celestial  bodies  at  the  beginning  of  the  twentieth  century,  the  science  of 
chaos  is  a  relatively  new  field.  It  is  one  which  essentially  became  possible  after 
the  birth  of  the  computer  which  not  only  allowed  numerical  integration  of 
systems  of  mathematical  equations,  but  also  shorter  and  shorter  time  intervals 
at  which  experimental  data  could  be  obtained.  In  the  1960's  Edward  Lorenz  of 
MIT  used  a  computer  to  numerically  integrate  twelve  differential  equations 
which  described  a  crude  model  of  the  weather.  In  his  simulations  he  found 
dramatic  differences  in  output  for  subtle  changes  in  input.  [Ref.  3:  p.  51]  Since 
then,  many  other  researchers  have  used  computers  to  investigate  chaotic 
behavior  in  virtually  every  field  including  mathematics,  biology,  chemistry, 
ecology,  economics,  physics,  medicine,  meteorology,  and  engineering.    James 


Gleick  presents  an  interesting  and  easy-to-read  development  of  the  science  of 
chaos  in  Chaos:   Making  a  New  Science  [Ref.  4]. 

The  chaotic  behavior  of  flexible  space  structures  has  been  studied  by 
Moon  and  Li.  In  their  analysis  of  a  3.5  meter  long  3-D  truss,  they  found  that  the 
free-play  in  the  joints  of  the  structure  resulted  in  chaotic  motion  and  modal 
frequencies  considerably  lower  than  those  determined .  by  linear  methods. 
Chaotic  dynamics  in  this  structure  and  other  space  structures  might  make  it 
difficult  to  design  active  controls  to  damp  out  transient  dynamics.  [Ref.  5]  In 
addition,  chaotic  dynamics  leads  to  inaccuracies  in  linear  finite  element  models 
which  are  currently  used  to  obtain  modal  analysis  during  the  design  of 
spacecraft.  The  results  of  the  Low  Power  Atmospheric  Compensation 
Experiment  (LACE)  may  lead  to  improvements  of  these  finite  element  models. 
[Ref.  6] 

In  another  study  by  Moon  and  Shaw,  chaotic  vibration  were  observed  in 
an  elastic  beam  with  non-linear  boundary  conditions.  For  certain  values  of 
forcing  frequency  and  amplitude,  the  periodic  motion  became  unstable  and 
nonperiodic.  The  results  were  typical  of  a  class  of  mechanical  oscillators  with 
play  or  amplitude  constraining  stops.  [Ref.  7] 

B.    GEOMETRIC  AND  TOPOLOGICAL  METHODS  OF  CHAOS 

The  purpose  of  this  section  is  to  describe  several  methods  which  have  been 
developed  for  detecting  chaotic  behavior  in  dynamical  systems.  The  methods 
which  are  discussed  are  time  series,  phase  plane,  frequency  domain,  Poincare 
section,  Lyapunov  exponent,  and  fractal  correlation  dimension.  The  first  four 
methods  are  more  qualitative  in  nature,  and  answer  the  question:  Is  chaos 
present?    The  last  two  methods,  Lyapunov  exponent  and  fractal  correlation 


dimension,  are  more  quantitative  in  nature.    Generally,  when  determining  if 
chaos  is  present  in  a  system,  several  of  the  above  methods  are  used  since 
applying  one  method  by  itself  is  not  always  indicative  of  chaotic  motion. 
1.    Time  Domain  Analysis 

The  traditional  way  of  observing  signals  is  to  view  them  in  the  time 
domain.  This  is  simply  a  record  of  what  is  happening  to  a  parameter  of  the 
system  versus  time.  Sometimes  the  first  indication  of  chaotic  motion  is  that  the 
parameter  observed  exhibits  no  pattern  or  periodicity.  This  test  is  not  always 
conclusive,  since  the  motion  observed  could  have  a  long  period  behavior  that  is 
not  readily  observable.  In  addition,  some  nonlinear  systems  may  appear  non 
periodic,  but  can  be  broken  down  into  several  periodic  signals.  [Ref.  2:  p.  42] 

Figure  2.1  illustrates  both  the  nonperiodic  nature  of  chaotic  signals  and 
also  the  sensitive  dependence  on  initial  conditions.  Thompson  and  Stewart 
numerically  integrated  the  Duffing's  equation 

x  +  0.05x  +  x3  =  7.5  cos  t  (2.1) 

to  produce  the  time  history  in  the  figure.  The  ragged  appearance  persists  for  as 
long  a  time  as  the  integrations  are  carried  out.  There  is  never  an  exact 
repetition  of  any  part  of  the  signal— it  is  truly  nonperiodic.  The  figure  also 
illustrates  that  for  slight  changes  in  initial  conditions,  the  trajectories  quickly 
diverge.  [Ref.  8:  p.  4] 


2.  Frequency   Domain   Analysis 

Because  of  the  problems  encountered  when  trying  to  detect  chaos  from 
the  time  domain,  the  frequency  domain,  or  Fourier  Spectrum,  is  primarily  used. 
Observing  the  signal  in  the  frequency  domain  will  instantly  reveal  if  the  signal  is 
made  up  of  a  sum  of  periodic  signals.  Chaotic  vibrations  will  exhibit  a  broad 
spectrum  of  frequencies  in  the  output  when  the  input  is  a  single  frequency 
harmonic  motion.  Subharmonics  and  Superharmonics  of  this  frequency  will 
appear  not  as  sharp  spikes  like  the  excitation  frequency  but  as  rounded  hills. 
Once  again,  this  test  is  not  always  conclusive  since  multiharmonic  outputs  may 
result  from  hidden  degrees  of  freedom  of  which  the  observer  is  unaware.  In 
order  to  detect  chaotic  vibrations  in  large  degree-of-freedom  systems,  it  is 
necessary  to  observe  changes  in  the  spectrum  while  one  parameter  such  as 
driving  amplitude  or  frequency  is  varied  [Ref.  2:  p.  47]. 

Figures  2.2  and  2.3  show  the  Fourier  spectrum  for  a  periodic  signal 
generated  by  the  equation 

x  =  3  cos  (47i  t)  (2.2) 

and  the  chaotic  signal  of  Equation  (2.1)  respectively.  Notice  that  Figure  2.2  has 
a  sharp  spike  only  at  the  excitation  frequency  while  Figure  2.3  has  multiple, 
broadband  spikes  at  frequencies  other  than  the  excitation  frequency  of  0.16  Hz. 

3.  Phase   Plane 

A  dynamical  system  consists  of  two  parts:  the  state  (essential 
information  about  a  system)  and  the  dynamic  (rule  that  describes  how  the  state 
evolves  with  time).   The  motion  of  the  system  can  be  visualized  in  a  state  space 


or  phase  plane.  The  phase  plane  is  a  construct  whose  coordinates  are 
components  of  the  state.  The  coordinates,  therefore,  depend  on  the  type  of 
system  being  studied.  A  mechanical  system  will  usually  use  position  versus 
velocity.  A  simple  example  to  describe  this  is  the  motion  of  a  pendulum.  All 
that  is  needed  to  determine  its  motion  are  position  and  velocity.  The  state  is 
thus  a  point  on  the  phase  plane.  Newton's  Laws  provide  the  dynamic, 
expressed  mathematically  as  a  differential  equation,  that  describes  how  the 
state  evolves.  Without  the  presence  of  friction,  the  pendulum  traces  out  a  circle 
or  orbit  in  the  phase  plane.  When  friction  is  present,  the  orbit  spirals  in  to  a 
single  point  or  attractor.  [Ref.  3:  p.  49]  This  describes  the  motion  of  a  simple 
non  chaotic  system.  The  orbits  of  chaotic  systems  never  close  or  retrace 
themselves.  Thus,  the  trajectory  of  the  orbits  will  tend  to  fill  up  a  section  of  the 
phase  plane.   Once  again,  this  is  not  always  a  conclusive  test  for  chaos. 

The  pseudo-phase  plane  is  another  technique  which  has  been  used 
when  only  one  variable  is  able  to  be  measured.  This  technique  was  discovered 
by  scientists  at  the  University  of  Santa  Cruz  while  analyzing  a  dripping  faucet. 
They  constructed  a  pseudo-phase  plane  using  a  measured  variable  and  a  time- 
shifted  version  of  the  same  measured  variable.  Takens  later  proved  that  a 
pseudo-phase  plane  exhibited  the  same  characteristics  as  a  phase  plane.  [Ref. 
4:  p.  266]  The  choice  of  the  time  shift  or  embedding  time  will  shrink  or 
expand  the  orbit  about  the  line  y  =  x.  Dvorak  and  Klaschka  have  proven  that 
an  embedding  time  between  five  and  eleven  time  intervals  produces  a  pseudo- 
phase  plane  with  the  least  amount  of  bias  [Ref.  9].  Chaotic  motions  exhibits  the 
same  characteristics  in  the  pseudo-phase  plane  as  they  do  in  the  phase  plane. 


The  pseudo-phase  plane  for  the  periodic  signal  appears  as  a  circle  as 
expected  (Figure  2.4).   This  signal  constantly  retraces  itself.    As  seen  in  Figure 
2.5,  the  pseudo-phase  plane  for  the  chaotic  signal  of  Equation  (2.1)  never 
retraces  itself  and  fills  up  a  section  of  the  plane. 
4.    Poincare   Section 

The  Poincare  section  is  constructed  by  passing,  a  two-dimensional 
surface  through  a  three-dimensional  phase  space  and  observing  where  the 
points  of  the  trajectory  intersect  the  plane.  This  slice  can  be  made  through  an 
infinite  number  of  locations  to  reveal  the  internal  structure  of  the  trajectory. 
Periodic  signals,  such  as  an  undamped  pendulum,  have  one  point  on  the 
Poincare  section  while  signals  with  multiple  harmonics  produce  a  finite  number 
of  points  on  the  Poincare  section.  If  the  Poincare  section  does  not  consist  of  a 
finite  set  of  points  or  a  closed  orbit,  the  motion  may  be  chaotic.  In  undamped  or 
slightly  damped  chaotic  systems,  the  Poincare  section  appears  as  a  cloud  of 
unorganized  points.  In  damped  chaotic  systems,  the  Poincare  section  appears 
as  an  infinite  set  of  highly  organized  points  which  look  like  parallel  lines.  If  after 
several  enlargements  the  structure  continues  to  exist,  the  system  is  described 
as  having  a  strange  attractor.  This  embedding  of  structure  within  structure  is 
referred  to  as  a  Cantor  set  after  the  mathematician  George  Cantor  (1845  - 
1918).  [Ref.  2:  p.  52] 

Figure  2.6  is  the  Poincare  section  of  the  periodic  signal.  Notice  that 
there  is  only  one  point  in  the  top  right  hand  corner.  The  damped,  chaotic  signal, 
however,  appears  as  a  set  of  points  forming  parrallel  lines  as  shown  in  Figure 
2.7.  If  this  figure  were  enlarged,  the  same  pattern  would  emerge  after  each 
enlargement.   This  is  referred  to  as  the  self-similar  property. 


5.    Lyapunov   Exponent 

The  Lyapunov  exponent  is  a  measure  of  the  average  rate  of 
divergence  of  trajectories  initially  separated  by  an  infinitesimal  amount.  A 
system  is  chaotic  if  it  has  one  or  more  positive  Lyapunov  exponents,  which 
means  nearby  trajectories  are  diverging  [Ref.  10:  p.  323].  Lyapunov  exponents 
have  proven  to  be  a  useful  dynamical  test  for  chaotic  systems.  Since  nearby 
trajectories  correspond  to  nearly  identical  states,  exponential  divergence 
between  the  trajectories  means  that  predictability  of  the  system  is  rapidly  lost. 
Systems  with  subtle,  unresolvable  initial  differences  quickly  behave  differently. 
[Ref.  1 1:  p.  285]  The  Lyapunov  exponent,  X,  is  defined  as 


,A.t 


or 


d(t)  =  d  2Kl  (2.3) 


^  =  }log2(^)  (2.4) 


where  :  d  is  the  initial  distance  between  the  two  trajectories 

d(t)  is  the  distance  at  a  later  time,  t 
As  indicated  in  Equation  (2.4)  a  positive  exponent  means  that    'd(t)'  is  greater 
than  'd'  and  the  trajectories  are  diverging.    A  zero  exponent  indicates  'd(t)'  and 
'd'  are  the  same  and  the  trajectories  coincide.    A  negative  exponent  indicates 
'd(t)'  is  less  than  'd'  and  the  trajectories  are  converging. 

A  system  with  n  degrees  of  freedom  has  a  set  of  n  exponents,  called 
the  Lyapunov  spectrum.  This  spectrum  characterizes  the  divergence  of  the 
trajectory  in  n  directions.  In  order  to  understand  the  meaning  of  the  Lyapunov 
spectrum,  imagine  a  sphere  of  neighboring  points  in  the  phase  space  about  a 
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central  point  on  the  trajectory.  As  the  points  evolve  with  time,  this  sphere  will 
deform.  If  one  axis  is  diverging  exponentially  the  sphere  stretches  along  that 
direction.  If  another  axis  is  converging,  the  sphere  shrinks  along  that  direction. 
Algorithms  are  available  to  calculate  the  exponents  of  the  spectrum  from  the 
measurement  of  a  single  observed  state  by  constructing  the  pseudo-phase 
space.  [Ref  12:  p.  89]  The  program  "Chaos"  uses  a  method  developed  by 
Wolfetal[Ref.  11]. 

6.    Fractal  Correlation  Dimension 

The  fractal  correlation  dimension  is  defined  as 

C(r)  =  rd  (2.5) 

where:  C(r)  is  the  probability  that  a  point  on  the  attractor  is  within  a 

circle,  sphere,  or  hypersphere  of  radius  V 

r  is  the  radius  of  the  circle,  sphere,  or  hypersphere 

d  is  the  fractal  correlation  dimension. 
Solving  for  'd'  in  Equation  (2.5)  yields 


,     ..     ,ln[C(r)], 
d  =  hm{lnfrr}  (2-6) 


C(r)  may  be  calculated  by  constructing  a  sphere  at  each  point  on  the 
trajectory  and  counting  the  number  of  points  in  each  sphere,  such  that 


C(r)  =  lim^-iiH(r-|xi-Xj|)  (2.7) 

N2    i     j 
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where:  H(s)  is  the  Heaviside  function  and  is  equal  to  one  if 's'  is  greater 

than  zero  and  equal  to  zero  if  Y  is  less  than  zero 
I  xi  -  xj  I  is  the  Euclidean  distance  between  points 
N  is  the  total  number  of  points  [Ref.  2:  p.217] 
The  procedure  used  in  the   "Chaos"  program  by  Sarigul-Klijn  to 

compute  the  fractal  correlation  dimension  is  described  below. 

1.  Start  with  a  point  on  the  attractor  and  calculate  the  number  of  points 
inside  a  circle  of  radius,  r. 

2.  Calculate  the  probability  C(r)  by  dividing  this  number  of  points  by  the 
total  number  of  points  in  the  attractor. 

3.  Repeat  this  for  several  points  along  the  attractor. 

4.  Compute  C(r)  for  several  values  of  r. 

5.  The  slope  of  log[C(r)]  versus  log(r)  yields  d,  the  fractal  correlation 
dimension,  for  the  chosen  number  of  pseudo-phase  space  dimensions. 
[Ref.  l:p.  88]. 

6.  Repeat  steps  1  through  5  for  higher  pseudo-phase  space  dimensions  to 
obtain  an  asymptotic  estimate  of  the  fractal  correlation  dimension  for  the 
dynamical  system. 

Nonchaotic  systems  have  integer  fractal  correlation  dimensions  such  as 
an  equilibrium  point  (zero),  limit  cycle  (one),  or  toroidal  attractor  (two).  The 
fractal  correlation  dimension  of  a  chaotic  system  is  non-integer  and  is 
independent  of  the  dimension  of  the  phase  space  used  for  the  calculation.  [Ref. 
13:  p.  101] 
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Figure  2.1    Divergence  from  adjacent   starts. 
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Figure  2.2   Fourier  spectrum  of  a  periodic  signal. 
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Figure  2.3   Fourier  spectrum  of  a  chaotic  signal. 
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PHASE    PLANE 
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Figure  2.4   Pseudo-phase  plane  of  a  periodic  signal. 
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Figure  2.6   Poincare  section  of  a  periodic  signal. 
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Figure  2.7   Poincare  section  of  a  chaotic  signal. 
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III.     THEORETICAL    ANALYSIS 

A.   DUFFING'S    EQUATION 

The  mathematical  equation  which  describes  the  motion  of  a  nonlinear 
elastic  beam  is  represented  by  the  well-known  Duffing's  equation. 

mv +  cv +  kiv  +  IC3V3  =  Acos  cot  (3.1) 

where:  m  =  mass  of  the  beam 

c  =  damping  coefficient 
ki  =  linear  stiffness  coefficient 
k3  =  nonlinear  stiffness  coefficient 
A  =  amplitude  of  the  excitation 
co  =  frequency  of  the  excitation 
v  =  transverse  displacement  of  the  beam 

Once  the  coefficients  of  the  Duffing's  equation  are  known  it  can  be  integrated 
numerically  to  solve  for  the  displacement,  velocity  and  acceleration  of  the  tip  of 
a  cantilever  beam.  A  fourth  order  Runge-Kutta  scheme  (see  Appendix  A)  was 
developed  to  perform  the  integration.  The  following  methods  are  used  to 
determine  the  coefficients  using  the  properties  of  a  beam  made  of  7075-76 
sheet  aluminum  listed  in  Table  3.1. 
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B.    ASSUMED   MODES   METHOD 

In  order  to  determine  the  coefficients  to  be  used  in  Duffings  equation,  the 
assumed  modes  method  is  used.  This  method  is  an  extension  of  the  principle 
of  virtual  displacements  and  is  used  to  produce  a  generalized  parameter  model 
that  approximates  the  flexible  behavior  of  the  beam.  A  virtual  displacement  is 
an  infinitesimal,  imaginary  change  in  the  displacement  of  a  system  as  shown  in 
Figure  3.1. 

The  dotted  curve  in  the  figure  shows  the  virtual  displacement,  5v(x,t),  of  the 
beam.  It  must  satisfy  the  same  boundary  conditions  as  the  actual  displacement, 
v(x,t)--the  solid  curve  in  Figure  3.1.  The  virtual  displacement  is  not  a  function 
of  time  in  the  same  sense  as  the  actual  displacement  but  is  rather  a  small 
change  of  configuration  relative  to  the  configuration  of  the  beam  at  the  time  "t" 
[Ref.  14:  p.  33]. 

To  create  the  single  degree  of  freedom  model,  the  deflection  of  the  beam 
may  be  approximated  by 

v(x,t)  =  \|<x)  fl(t)  (32) 

Similarly,  the  virtual  displacement  is  given  by 

6v(x,t)  =  \j/(x)  5d(t)  (3.3) 

The  shape  function,  \j/(x),  may  be  any  admissible  function  which  adequately 
describes  the  behavior  of  the  system  and  possesses  derivatives  of  an  order  at 
least  equal  to  those  specified  by  the  boundary  conditions.      The  generalized 
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displacement,  d(t),  is  the  solution  to  Duffing's  equation.  [Ref.  14:  p.  34]  The 
shape  function  which  is  chosen  is  the  characteristic  function  for  a  clamped-free 
beam. 

V(5)  =  VJ  sin  <«5  -*)  +  exp  (-J5)  +  exp  ( - 1  (1  -  £))  (3#4) 

where  ^  is  the  beam  coordinate,  nondimensionalized  with  respect  to  the  length, 
L,  and  varying  between  0  and  1.  [Ref  15:  p.  439] 

The  boundary  conditions  for  a  cantilever  beam  are 


BC1 
BC2 
BC3 
BC4 


v(0,t)  =  0 
v'(0,t)  =  0 
v"(U)  =  0 
v'"(L,t)  =  0  V-y) 


where  (')  indicates  the  derivative  with  respect  to  the  spatial  coordinate  x. 

The  first  boundary  condition  states  that  there  is  no  displacement  at  the  fixed 
end  while  the  second  states  that  its  slope  must  also  vanish  here.  The  third  and 
fourth  boundary  conditions  state  that  there  must  be  no  moment  and  no  shear 
forces  at  the  free  end.  Applying  the  boundary  conditions  of  Equation  (3.5)  to 
the  shape  function  of  Equation  (3.4)  gives  only  a  0.7%  error  in  amplitude  at  the 
midpoint.  [Ref  16:  p.  1014] 

The  first  and  second  derivatives  of  the  shape  function  will  be  necessary  in 
determining  the  coefficients  for  Duffing's  equation.   Subsequently, 


V^)=V2|cos(^-^)-S-exp(-^)  +  2.exp(-|(l-^))  (3>6) 

¥"©  =  -VI^-sin(^-^)+^-exp(-^)  +  ^-exp(-|(l-5))        (3.7) 
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The  principle  of  virtual  displacements  states  that  the  combined  virtual  work 
of  all  forces  must  be  zero.  More  explicitly  this  is  stated  as 

6Wnc  -  5V  +5Wdamp  +  SW^  =  0  (3.8) 

where  8Wnc  is  the  virtual  work  of  nonconservative  forces,  8V  is  the  change  in 
strain  energy  of  the  system,     5Wdamp    is  the  virtual  work  of  damping  forces, 
and  5Winert  is  the  virtual  work  of  inertia  forces  [Ref.  14:  p.  34]. 
1.    Mass   Coefficient   Determination 

Solving  for  the  virtual  work  of  inertia  forces  will  yield  the  mass  term,  m, 
of  Duffing's  Equation  (3.1).  Since  the  inertia  force  per  unit  length  is  distributed 
along  the  beam,  the  virtual  work  of  inertia  forces  is 

5Winert  =1     "PA  V  5v  dx 

JO 

=  -pAL  J   \^d^5d  (3.9) 

Jo 

Using  the  values  for  density,  area,  and  length  of  the  beam  contained  in  Table 
3.1  and  numerically  integrating  using  MATHCAD  (see  Appendix  B),  the  work 
done  by  inertia  forces  is 

SW^t  =  1.0360  x  101  £  8#  (3.10) 

Therefore,  the  mass  coefficient,  m,  of  Duffing's  equation  is  0.10360  slugs. 
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2.    Stiffness   Coefficients   Determination 

If  a  beam  is  modeled  as  shown  in  Figure  3.2  with  coordinate  axes  x,  y, 
and  z  and  corresponding  displacements  u,  v,  and  w,  it  can  be  shown  [Ref.  17:  p. 
22]  that  the  strain  in  the  x  direction  is  given  by  the  following  equation: 


_  du     i 
dx      2 


dx  dx  dx 


(3.11) 


Assuming  that  all  higher  order  terms  are  negligible  except  for  (dv/3x)2  (small 
strain,  large  rotation  assumption),  the  above  equation  becomes: 


du      i  ,9v,2 


Applying  the  classical  Euler-Bernoulli  assumption  that  cross  sections  are 
assumed  to  remain  planar  and  normal  to  the  centroidal  axis  after  deformation 
[Ref.  18:  p.  152],  it  can  be  seen  from  Figure  3.3  that 


u  =  -ytan0z  (3.13) 

_Av 


tan  0Z  = 


Ax 


(3.14) 


Subsequently 


u=-y^  0.15) 

du  ci2v 

§r-y^  (3-i6) 
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Substituting  Equation  (3.16)  into  Equation  (3.12)  yields 


exx  =  -y^  +k—)  (3  17) 

dx2      2dx2'  KD      } 


The  expression   for   the   strain   energy   of  a  beam   undergoing   transverse 
deflection  is  given  by  [Ref.  19:  p.  169] 


V  =  ±l    cxxexxdxdA  <3-18) 


Jo 


Because  the  beam  under  study  is  not  deflected  so  the  material  passes  its  yield 
point,  it  behaves  as  a  Hookean  material  such  that 

oxx  =  Eexx  (3.19) 

Therefore  the  strain  energy  becomes 

V=1E|   e?,dxdA  (3-2°) 


-i-f 

JO 


Squaring  Equation  (3.17)  and  substituting  it  into  Equation  (3.20)  yields: 


fy2(d^v)2  .  y(&L)2&!L  +  J^dv)^ 
[y  \ix2;      yW  dx2     4\lx' 

o 


dxdA 


(3.21) 
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If  the  intersection  of  the  y  and  z  axes  lies  at  the  centroid  then 


'. 


y2  dA  =  Izz  (3.22) 

ydA  =  0  (3.23) 


Therefore,  Equation  (3.21)  becomes: 


v-¥f&>2**¥JV*  (124) 

Jo  Jo 


Applying  variational  calculus,  Equation  (3.24)  becomes 


5V  =  EI  f  dWSv        +EAf  (dv}3d5v  ^  (3.25) 

I  dx2  dx2  2  j    W     dx 


which  can  be  written 


8V  =  E  I  J   v"  (5v)"  dx  +  EA  J    (V')3(5v)'  dx  (3-26) 

Jo  Jo 

Substituting  the  derivatives  of  Equation  (3.2)  and  (3.3)  and  substituting 
Equations  (3.6)  and  (3.7)  into  Equation  (3.26)  yields 

5V  =  EI      (\|T)2  ■&  Sfl  dx  +  EA-  I    (y')4  #3  8#  dx  ^21"> 

Jo  Jo 


25 


Substituting  in  the  values  from  Table  3.1  and  numerically  integrating  using 
MATHCAD  (See  Appendix  B),  Equation  (3.27)  becomes 

5V  =  1.5433  x  10"2  fl  8fl  +  4.5943  x  101  fl3  5d  (3.28) 

From  Equation  (3.28),  the  linear  and  nonlinear  stiffness  coefficients,  kiand  k3, 
are  seen  to  be  0.015433  lbs/in  and  45.943  lbs/in^  respectively. 
3.    Damping  Coefficient   Determination 

If  the  beam  is  assumed  to  be  modelled  by  a  linear  single  degree  of 
freedom  system  the  damping  coefficient  is  given  by 

cd  =  2Vk1m  £  (3.29) 

where:  £  =  viscous  damping  factor 

The  values  for  ki  and  m  were  determined  previously.  In  order  to 
determine  the  viscous  damping  factor,  the  logarithmic  decrement  method  was 
used  [Ref.  14:  p.  62].  The  logarithmic  decrement,  5,  is  defined  by 


6.  to  #).**££  (3.30) 


where:        Up  =  amplitude  of  motion  at  beginning  of  cycle 
Uq  =  amplitude  of  motion  after  N  cycles 
N  =  number  of  cycles 
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The  above  equation  is  used  to  solve  for  the  viscous  damping  factor  by 
displacing  the  tip  of  the  beam  and  then  measuring  the  two  amplitudes,  Up  and 
Uq,  and  counting  the  number  of  cycles  between  the  two  measurements.  These 
measurements  were  performed  on  an  air  table  in  the  Flexible  Spacecraft 
Simulator  Laboratory  of  NPS  as  shown  in  Figures  3.4  and  3.5.  The  beam  was 
horizontally  mounted  on  the  table  and  an  air  pad  was  attached  to  the  tip  of  the 
beam.  The  tip  was  then  deflected  and  the  resulting  motion  achieved  essentially 
without  friction.  The  values  for  Up,  Uq,  and  N  appear  in  Table  3.2.  Using 
Equation  (3.30),  the  average  viscous  damping  factor,  £,  was  determined  to  be 
6.67x10*3  or  approximately  0.7%. 

The  virtual  work  of  the  damping  forces  is  given  by 


( 


5Wdamp  =  -  I    cd  v  5v 

=  -cd  I    ¥2d^-3  5i3  (3.31) 


Substituting  in  the  value  for  cd  obtained  from  Equation  (3.29)  and  numerically 
integrating  once  again  yields 

5Wdamp  =  2.5654  x  10"3  0  5i5  (3.32) 

Therefore,  the  damping  coefficient,  c,  of  Equation  (3.1)  is  0.0025654  slugs/s. 
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C.    THEORETICAL   RESULTS 

Appendix  A  contains  the  FORTRAN  computer  code  used  to  numerically 
integrate  the  Duffing's  equation  which  mathematically  describes  the  motion  of 
the  flexible  beam  being  studied.  Table  3.3  summarizes  the  results  of  the 
previous  section  in  which  the  coefficients  of  Equation  (3.1)  were  determined. 
Equation  (3.1)  was  then  integrated  using  an  amplitude  of  3  lbs.  (the  force 
applied  by  the  shaker)  and  a  frequency  of  5.6  Hz.  The  initial  displacement  and 
velocity  chosen  were  1  inch  and  0  in/s,  respectively.  The  results  were  then 
formatted  to  be  read  by  the  computer  program  "Chaos"  developed  by  M. 
Sarigul  Klijn  [Ref.  1]. 

Qualitative  tests  for  chaos  are  satisfied  by  the  time  series,  Fourier  spectrum, 
pseudo-phase  plane  and  Poincare  section.  Figure  3.6  shows  the  time  history  of 
the  acceleration  of  Equation  (3.1).  It  shows  no  noticeable  periodic  behavior  and 
this  is  confirmed  in  the  Fourier  spectrum  of  Figure  3.7.  The  power  spectrum  is 
made  up  of  many  frequencies  with  broad  band  characteristics.  The  pseudo- 
phase  plane  of  Figure  3.8  was  obtained  using  an  embedding  time  of  10.  This 
embedding  time,  consistent  with  the  findings  of  Dvorak  and  Klaschka  [Ref  9], 
yielded  the  least  amount  of  bias.  Notice  that  the  trajectory  frequently  crosses 
itself  and  fills  up  a  section  of  the  phase  space.  The  Poincare  section,  as  shown 
in  Figure  3.9,  is  also  characteristic  of  damped  chaotic  systems  with  its  cloud  of 
unorganized  points  filling  a  section  of  the  space.  It  was  obtained  by  passing  a 
plane  through  the  three-dimensional  pseudo-phase  space  at  z  =  0. 

The  quantitative  tests  for  chaos  were  also  performed.  Figure  3.10  is  a  plot 
of  the  average  Lyapunov  exponent.  This  figure  is  obtained  by  time  averaging 
the  many  computations  of  Equation  (2.4)  and  then  plotting  them  versus  the 
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number  of  samples,  4096.    The  average  Lyapunov  exponent  for  this  system  is 
0.949.   This  exponent  is  indicative  of  chaotic  systems. 

Figure  3.11  shows  the  plot  from  which  the  fractal  correlation  dimension  is 
determined.  As  mentioned  previously,  the  fractal  correlation  dimension  is  the 
asymptotic  value  obtained  by  plotting  the  slope  of  log[C(r)]  versus  log(r)  for 
successively  higher  pseudo-phase  space  dimensions.  It  was  determined  to  be 
approximately  2.4.  Moon  reports  a  value  of  2.5  for  his  buckled  beam 
experiment.  [Ref.  2,  p.  232].  The  non-integer  fractal  correlation  dimension  is 
yet  another  indication  of  chaos  in  dynamic  systems. 
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TABLE  3.1 


Beam   Properties 


Density,  p 

0.101  lbm/in3 

Modulus,  E 

10.4xl06psi 

Length,  L 

90.0  in 

Area,  A 

0.24  in2 

Thickness,  t 

0.080  in 

Mom  Inert,  I 

0.000128  in4 

Width,  w 

3.0  in 

TABLE  3.2            Logarithmic   Decrement   Values 

N 

Up  (in) 

U0  (in) 

c 

20 

2.5625 

1.0000 

0.00749 

20 

2.3125 

0.9375 

0.00718 

20 

2.2500 

0.9375 

0.00697 

20 

2.6250 

1.0625 

0.00720 

20 

2.7500 

1.1250 

0.00711 

25 

2.2500 

0.8125 

0.00648 

25 

2.6250 

0.8750 

0.00699 

25 

2.5000 

0.8750 

0.00668 

25 

2.2500 

0.8125 

0.00648 

25 

2.2500 

0.8125 

0.00648 

30 

1.8125 

0.6250 

0.00565 

30 

2.6250 

0.7500 

0.00665 

30 

2.1875 

0.6875 

0.00614 

30 

1.9375 

0.5625 

0.00656 

30 

1.9375 

0.6250 

0.00600 

30 


TABLE  3.3  Coefficients  of  Duffing's  Equation 


m  = 

=  1.0360  x 

io-i 

c  = 

:  2.5654  x 

10-3 

kl  = 

=  1.5433  x 

10-2 

k3 

=  4.5943  x 
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Figure  3.1    Illustration  of  virtual  displacements 
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Figure  3.2   Coordinates  of  a  simple  beam 
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Figure  3.3   Kinematics  of  a  simple  beam 
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Figure  3.4   View  1  of  beam  on  air  table 
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Figure  3.5   View  2  of  beam  on  air  table 
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Figure  3.6   Time  series  of  Duffing's  Figure  3.7    Fourier  spectrum 

Equation.  Duffing's  equation. 
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Figure  3.8   Pseudo-phase   plane   of 
Duffing's  Equation. 
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Figure  3.9   Poincare    section 
of  Duffing's  equation. 
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Figure  3.10  Lyapunov  Exponent  of  Duffing's  equation. 
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Figure  3.11 


Fractal  Correlation  Dimension  of  Duffing's  equation 
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IV.     EXPERIMENTAL    ANALYSIS 

A.  EXPERIMENTAL    SETUP 

The  aluminum  beam  with  properties  listed  in  Table  3.1  is  instrumented  with 
strain  gauges  and  an  accelerometer  as  shown  in  Figure  4.1.  The  type  of  strain 
gauges  used  are  one-quarter  inch,  350  ohm,  straight  geometry  gauges  from  the 
Micro  Measurements  Group.  The  accelerometer  is  a  Sensotec  model  JTF  flat 
pack.  This  accelerometer  has  a  weight  of  one  ounce,  a  peak  range  of  ±5  g's,  a 
sensitivity  of  10  mV/g,  and  a  usable  frequency  range  of  DC  to  150  Hz. 

The  beam  is  mounted  vertically  to  an  MB  Dynamics  Model  PM500A 
Vibration  Exciter.  The  outputs  from  the  strain  gauges  and  accelerometer  are 
then  fed  into  a  bridge  balancing  box  and  amplifier.  The  bridge  box  is  capable  of 
handling  four  inputs,  one  from  the  accelerometer  and  three  from  the  strain 
gauges.  From  the  bridge  box,  the  signal  goes  to  an  analog-to-digital  converter. 
From  there  the  acceleration  and  strain  data  is  stored  on  a  Zenith  computer. 
Mechanical  stops  placed  58  inches  up  from  the  base  prevent  the  beam  from 
buckling  past  one  and  three-fourths  inches  on  either  side.  (Figure  4.2) 

B.  EXPERIMENTAL   PROCEDURE 

The  accelerometer  is  placed  at  location  Al  in  Figure  4.1.  The 
accelerometer  and  the  three  strain  gauges  are  then  connected  to  the  data 
acquisition  system.  The  beam  is  held  in  the  upright  vertical  position  in  order  to 
zero  out  all  the  readings  from  the  gauges  and  the  accelerometer  using  the 
bridge  balancing  circuit.  The  beam  is  then  released  while  simultaneously 
starting  the  exciter.   The  exciter  provides  an  excitation  frequency  of  5.6  Hz  and 
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a  displacement  amplitude  of  one-quarter  inch  peak-to-peak.  This  procedure  is 
then  repeated  twice,  placing  the  accelerometer  at  locations  A2  and  A3.  Strain 
gauge  measurements  at  location  SI,  S2,  and  S3  are  obtained  for  each  run. 

C.   EXPERIMENTAL    RESULTS 

The  data  from  each  of  the  runs  is  transferred  into  a  format  readable  by  the 
"Chaos"  program.  It  is  then  analyzed  using  the  techniques  discussed  in 
Chapter  II.  These  results  appear  at  the  end  of  this  chapter.  As  shown  in 
Figures  4.3  through  4.10,  the  four  qualitative  tests  for  chaos  seem  to  be  met: 

1.  Erratic  and  unpredictable  time  histories  from  both  accelerometers  and 
strain  gauges. 

2.  Multiple  broadbanded  power  spectrum  spikes. 

3.  Crossing  trajectories  and  filled  up  pseudo-phase  space. 

4.  A  cloud  of  unorganized  points  in  the  Poincare  section  characteristic  of 
lightly  damped  systems. 

These  results  appear  similar  to  the  chaotic  results  obtained  by  Sarigul-Klijn 
[Ref.  1]  in  his  analysis  of  helicopter  flight  test  data  using  higher  harmonic 
control. 

The  quantitative  tests  show  that  the  location  of  the  type  of  sensor  chosen 
impacts  whether  or  not  chaotic  motion  is  found.  Strain  gauges  only  indicate 
chaos  if  placed  near  the  root,  while  accelerometers  only  indicate  chaos  if 
placed  near  the  tip.  Figures  4.11  through  4.14  are  plots  of  the  average 
Lyapunov  exponents  and  fractal  correlation  dimensions.  A  summary  of  the 
Lyapunov  exponents  and  fractal  correlation  dimensions  for  each  sensor  appear 
in  Tables  4.1  and  4.2.  The  Lyapunov  exponents  are  all  positive  for  the  strain 
gauges.  For  the  accelerometers,  the  Lyapunov  exponent  diminishes  as  the 
accelerometer  is  placed  closer  to  the  root.     This  indicates  that  measuring 
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accelerations  near  the  root  reduces  the  chances  of  obtaining  chaotic  results. 
This  is  anticipated  since  the  beam  does  not  experience  large  deflections  in  this 
area  and  can  be  characterized  by  linear  equations.  Likewise,  the  fractal 
correlation  dimensions  of  the  sensors  indicate  that  accelerometers  should  be 
placed  closer  to  the  tip  and  strain  gauges  closer  to  the  root.  The  accelerometer 
nearest  the  root  has  an  integer  fractal  correlation  dimension,  like  that  of  a  non- 
chaotic  system.  The  two  nearer  to  the  tip  have  non-integer  fractal  dimensions, 
characteristic  of  chaotic  systems.  Conversely,  the  two  strain  gauges  near  the 
root  have  non-integer  fractal  dimensions,  while  the  one  closer  to  the  tip  has  an 
integer  fractal  dimension. 
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TABLE  4.1           Lyapunov   exponents 

SI 

0.272 

Al 

0.364 

S2 

0.237 

A2 

0.0397 

S3 

0.390 

A3 

-0.0282 

TABLE  4.2 


Fractal  Correlation  Dimensions 


SI 

4.8 

Al 

4.4 

S2 

4.4 

A2 

4.4 

S3 

4.0 

A3 

4.0 
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Figure  4.2   Experimental   Set-up 
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Figure  4.3   Time  series  of  strain  gauges 

(a)  SI      (b)  S2     (c)  S3 


42 


TIME  SERIES 


TME  SERIES 


(c) 
Figure  4.4   Time   series   of  accelerometers 
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Figure  4.5   Pseudo-phase  planes  of  strain  gauges 
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(c) 
Figure  4.6   Pseudo-phase   planes   of  accelerometers 
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Figure  4.7   Fourier  spectrums  of  strain  gauges 
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Figure  4.8   Fourier   spectrums   of  accelerometers 
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Figure  4.10 
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V.   CONCLUSIONS  AND  SCOPE  FOR  FUTURE 

RESEARCH 

For  the  first  time  geometric  and  topological  methods  of  chaos  were  applied 
to  data  obtained  simultaneously  by  accelerometers  and  strain  gauges.  It  was 
determined  that  data  from  both  accelerometers  and  strain  gauges  provided 
similar  indications  to  the  presence  of  chaos  when  analyzed  using  time  series, 
Fourier  spectrums,  pseudo-phase  planes,  and  Poincare  sections.  Using  these 
techniques,  both  types  of  sensors  indicated  chaos  without  regard  to  the  location 
of  the  sensor.  However,  the  Lyapunov  exponents  and  fractal  correlation 
dimensions  indicated  that  chaotic  motion  is  only  detected  by  accelerometers  if 
they  are  placed  near  the  tip  and  that  chaotic  motion  is  only  detected  by  strain 
gauges  if  the  are  placed  near  the  root. 

The  presence  of  chaos  in  highly  flexible  spacecraft  structures  was  also 
investigated.  It  was  determined  that  the  geometric  nonlinearities  of  these 
structures  may  result  in  chaotic  vibrations  and  dynamic  instabilities.  These 
instabilities  may  then  affect  control  mechanisms  of  the  spacecraft.  Further 
analysis  could  be  performed  in  using  the  geometric  and  topological  methods  of 
chaos  to  analyze  the  effectiveness  of  these  mechanisms  to  adequately  control 
the  spacecraft  and  provide  improvements  into  the  way  the  control  is  performed. 
Chaotic  methods  may  also  have  an  application  in  updating  the  linear  finite 
elements  models  of  flexible  space  structures  with  more  accurate  nonlinear 
models. 
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APPENDIX 


TXTLSl 

DurriMO 

BATEi 

07  DEC  1001 

AUTHOR: 

EOBEBT  O.  VAUQBAJf 

STBTEXt 

VAX 

* 
* 
* 

*  COUP I LEI:       MAS* 

*  DESCXXPTXOMl 

*  This  progran  selves  Duff in?' ■  equation 
t 


*  ■*<T2/dt*2U)  ♦  c*d/dtU>  ♦  kl*s  ♦  k3*x'3  -  Acos(w*t) 

*  by  transforming  it  into  tvo  ordinary  differential  equations 

* 

*  d/dt{«)  -  y 

*  d/at(y)  «  Acos(v«t)  -  c«y  -  «1»«  -  *3«x"J 

* 

*  and  using  a  fourth  order  Stunga— Kufcta  SohaM  to  integrate  to 

*  find  ths  displacement,  x.   and  tba  velocity,  y. 

* 

*  The  amplitude  and  frequency  of  tha  forcing  funotion,  tha  etep 

*  aiia  for  intsgration,  tha  number  of  integration  steps,  and  the 

*  initial  values  foe  displacement  and  velocity  are  entered  in 

*  that  order  In  a  data  file  titled  "DUrrXMlT."  The  data  file 

*  titled  'Durrour"  contains  the  integrated  values  for  displacement 

*  and  velocity  in  a  format  compatible  vlth  the  program  "CBAOS" 

*  by  H.  Barigul-Elijn. 

* 

*  To  alter  this  program  to  eolve  an  aquation  other  than  Duffing' s 

*  above,  the  parameter,  HI,  must  be  changed  to  correspond  to  the 

*  number  of  ODB's  end  the  subroutine  "DiPfEQ"  sunt  be  changed 

*  to  include  the  proper  equation*. 

a 

*  Two  subroutines  ara  eaployad. 

*  DirriQ  calculates  the  values  for  the  ODE's 

*  RK4  is  the  Jtunge-Kutte  routine 

»  As  written  the  program  does  not  use  an  adaptive  step-site 

a 

*  The  var labia*  are  defined  aa  follow* i 

*  HE  -  the  number  of  ODE's  in  the  system 

*  MS  -  tha  number  of  steps  (iterations)  to  be  performed 

*  B  -  the  desired  etepsise 

*  A  -  amplitude  of  the  forcing  function 

*  F  -  frequency  of  the  forcing  function 

*  TIME,  T  -  time 

*  I,  XNDBX  -  integer  counters 

*  X  -  an  array  corresponding  to  the  independent  variables, 

*  I.e.  *(l)-x,  x(2)-t,  etc. 

*  xdot  -  an  array  corresponding  to  the  ODB'a 


PKOGRAM  DUrriMG 

PAAAMSTEB  (MB-2) 

IMPLICIT  DOUBLE  PBECI£IOH(A-H,X,M,0-Z) 

CKAXACTEB"7    DuTUTT 

DIMENSION   X(NE),XS0T(NB) 
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0PEN(7,riLE-'DUrriNIT' , STATUS- 'OLD') 

open(8,fxle-'duffout* , status-' old' ) 

TIME-0.0D0 
DUMMY-'  A12  3  4A' 

RBAD(7,*)  h,r,H,N5 
DO  5  I-1,KE 

READ(7,*)  X(I) 

5  CONTINUE 

WRITE(8,6)  A,T 
WRITR(8,7)  H,NB 
WRITE(8.8)  X(1).X(2) 

6  FORMAT  (IX, 'FORCING  FUNCT  AMPLITUDE  -  ' , F4 . 1 , 5X, ' FREQ  -  ' , F6 . 3 ) 

7  FORMAT  (IX,* STEP  3I2E  -  '  , r6 . 4 , 5X,  ' NUKBER  Or  STEPS  -  ',15) 

8  FORMAT  (1X,'INIT  DISP  -  ' , P4 . 1 , 5X, ' IKIT  VRL  -  »,F4.1) 

CALL  DIFrEQ(X,XDOT,A,F,TIME) 

WRITE( 8, 100)  DUMMY,TIME,X(1),X(2) ,XDOT(2) 

DO  15  I-1,N3 

DO  10  INDEX-=1»4 

CALL  DIFFEQ(X,XDOT,A,F,TIME) 
CALL  RK4 ( TIME, X,XDOT,H, INDEX) 

10  CONTINUE 

WRITE(8,100)  DUMMY,TIME,X(l),X(2),XDOT(2) 

15     CONTINUE 

100     FORMAT(A7,  F6.3,  Ell. 4,  2E12.4) 

STOP 

END 

************************************************************************ 

SUBROUTINE  RK4(T,X,XDOT,H, INDEX) 

*  This  subroutine  performs  numerical  integration  using  4th  order 

*  Runge-Kutta  methodc.   It  was  developed  by  I.M.  Robe. 

PARAMETER (NE- 2) 

IMPLICIT  DOUBLE  PRECISION! A-H, K,M,0-Z ) 

DIMENSION  X(NE) ,XDOT(NE) ,8AVED(NE) ,SAVEX(NE) 

GO  TO  (1,2,3.4) .INDEX 

1     DO  10  I-1,NE 

3AVEX( I )-X( I ) 

SAVED(I)-XDOT(I) 
10     X(I)-SAVEX(I)+0.5D0*H*XDOT(I) 
T-T+0.5D0*H 
RETURN 
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2  DO  20  I-1,NE 

SAVED ( I ) -SAVED ( I )  +  2  .  ODO *XDCT ( I ) 
20     X(I )«SAVEX(I )+0.5D0*H*XDOT(I) 
RETURN 

3  DO  30  I-1,KE 

SAVED(I)-8AVED( I )+2 . 0D0*XDOT( I ) 

30      X< I )-SAVEX( I )+H*XDOT( I ) 

T-T+0.5D0*H 
RETURN 

4  DO  40  I-1,NE 

40     X(I)-SAVBX(I)4-H/6.0D0*(SAVKD(I)+XDOT(I)  ) 
RETURN 

END 

********************** *************4*4^*** ************************* ****** 

SUBROUTINE  DIFFEQ(X,XDOT,A, F,T) 

*  This  subroutine  obtains  values  for  the  ode's  for  Duffing' 6 

*  equation  and  oust  be  changed  for  a  new  system  of  equations. 

PARAMETER!'  NE""2  ) 

IMPLICIT  DOUBLE  PRECISION( A-H,K,M,0-Z ) 

DIMENSION  X(NE) ,XD0T(NE) 

PI-4 . 0D0*DATAN{ 1 . ODO ) 

*  Be  Careful  of  units!    Kl  -  LBf/in.    K3  -  LBf/in*3 

*  M  »  slugs       C  -  slugs/s 
M-1.0360D-01 

C-2.5654D-03 

K1-1.5433D-02 

K3-4.5943D+01 

*  Divide  M  by  12  in/ft  so  that  XD0T(2)  has  units  in/s"2 

M-M/12.0D0 
XDOT(l)-X(2) 

XDOT(2)-(A*COS(F*2.0D0*PI*T)-C/12.0D0*X(2)-Kl*X(l)-K3*X(l)**3)/M 

RETURN 
END 
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APPENDIX         B 
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